# -*- coding: utf-8 -*-
"""
梯度上升发解决主成分分析问题
Created on Tue Feb 27 19:53:07 2018

@author: Allen
"""

import numpy as np
import matplotlib.pyplot as plt

X = np.empty( ( 100, 2 ) )
X[ :, 0 ] = np.random.uniform( 0., 100., size = 100 )
X[ :, 1 ] = 0.75 * X[ :, 0 ] + 3. + np.random.normal( 0., 10., size = 100 )

plt.scatter( X[ :, 0 ], X[ :, 1 ]  )
plt.show()


def demean( X ):
    return X - np.mean( X, axis = 0 )

X_demean = demean( X )

plt.scatter( X_demean[ :, 0 ], X_demean[ :, 1 ]  )
plt.show()

#print( np.mean( X_demean[:, 0] ), np.mean( X_demean[:, 1] ) )

# 梯度上升法
def f( w, X ):
    return np.sum( ( X.dot( w ) ** 2 ) ) / len( X )

def df_math( w, X ):
    '''求梯度'''
    return  X.T.dot( X.dot( w ) ) * 2. / len( X )

def direction( w ):
    res = w / np.linalg.norm( w )
    return res

def gradient_ascent( df, X, initial_w, eta, n_iters = 1e4, epsilon = 1e-8 ):
    '''梯度上升法'''
    w = direction( initial_w )# 因为求的是方向向量，所以需要单位化
    cur_iter = 0
    while cur_iter < n_iters:
        gradient = df( w, X )
        last_w = w
        w = w + eta * gradient
        w = direction( w ) # 注意：每一次求出w，都要向量单位化
        
        if abs( f( w, X ) - f( last_w, X ) ) < epsilon:
            break
        
        cur_iter += 1
    return w

# 初始化 w， 注意：不能使用 0，不然 函数就永远为0，所以需要随机出来一个
initial_w = np.random.random( X.shape[1] )

eta = 0.001

# 注意： 不能使用 StandardScaler标准化数据(归一化)
# 原因在于，pca本身就是求的方差最大化，如果先归一化了，方差也就确定了，pca就没意义了

w = gradient_ascent( df_math, X_demean, initial_w, eta )
print( w )

# 画出单位向量的图
plt.scatter( X_demean[:, 0], X_demean[:, 1] )
plt.plot( ( 0, w[0] * 50 ), ( 0, w[1] * 50 ), color = "r" )
plt.show()

# 查看一个极端情况
X2 = np.empty( ( 100, 2 ) )
X2[ :, 0 ] = np.random.uniform( 0., 100., size = 100 )
X2[ :, 1 ] = 0.75 * X2[ :, 0 ] + 3.
plt.scatter( X2[ :, 0 ], X2[ :, 1 ]  )
plt.show()

X2_demean = demean( X2 )
initial_w = np.random.random( X.shape[1] )
eta = 0.001
w1 = gradient_ascent( df_math, X2_demean, initial_w, eta )
plt.scatter( X2_demean[:, 0], X2_demean[:, 1] )
plt.plot( ( 0, w1[0] * 50 ), ( 0, w1[1] * 50 ), color = "r" )
plt.show()

'''
本代码中心思想：
    求出样本第一个主成分
    将两维度映射到一维，使得各个样本点之间方差最大
'''